
(NASA-CR-3990) TURBULENT E10W HELD 
PREDICTIONS IN SHARPLY CURVED IUEN-AROUND 
DUCTS Final Report (NASA) 34 p 

HC AO 3/M F A0 1 CSCL 01A Unclas 

H 1/0 2 43 159_ 


/ xmm 




| |j r | , p - pj 





NASA 



NASA Contractor Report 3990 


Turbulent Flow Field 
Predictions in Sharply 
Curved Turn-Around Ducts 


L. Michael Santi 

Memphis State University 
Memphis , Tennessee 


Prepared for 

George C. Marshall Space Flight Center 
under Grant NGT-0 1-008-021 


f\l/\SA 

National Aeronautics 
and Space Administration 

Scientific and Technical 
Information Branch 


1986 


TABLE OF CONTENTS 


Page 


INTRODUCTION 1 

OBJECTIVES 3 

BACKGROUND INFORMATION 5 

RESULTS 8 

CONCLUSIONS 22 

REFERENCES 25 


PRECEDING PAGE BLANK 


LIST OF ILLUSTRATIONS 


Figure Title Page 

1. Early SSME fuel side turbopump TAD design... 2 

2. Simple TAD analysis geometry 4 

3. Geometry and grid pattern for BFC TAD analyses: 20 x 80 grid 4 

4. Streamwise velocity profiles for simple TAD geometry 

with plug flow inlet: Analysis 1 10 

5. Pressure profiles for simple TAD geometry with plug 

flow inlet: Analysis 1 10 

6. Turbulent kinetic energy profiles for simple TAD geometry 

with plug flow inlet: Analysis 1 11 

7. Streamwise velocity profiles for simple TAD geometry with 

FDP flow inlet: Analysis 2 11 

8. Pressure profiles for simple TAD geometry with FDP flow 

inlet: Analysis 2 13 

9. Turbulent kinetic energy profiles for simple TAD geometry 

with FDP flow inlet: Analysis 2 13 

10. Streamwise velocity profile comparisons for simple 

TAD geometry 14 

11. Pressure profile comparisons for simple TAD geometry 14 

12. Pourahmadi and Humphrey model [9,10] staggered 

grid basis cell 15 

13. Behavior of Ff^C^) polynomial at various cross stream 

locations for 0 = 168 deg: Analysis 3 15 

14. Characteristic profiles for simple TAD geometry with plug 

flow inlet condition using the modified k-e model: Analysis 3 17 

15. Streamwise variation of enclose to the concave wall: Analysis 3 17 

16. Turbulent kinetic energy profiles for simple TAD geometry 

using the modified k-e model: Analysis 3 18 

17. Comparison of near wall turbulent kinetic energy predictions 18 

18. TAD I/O velocity profiles for BFC computation: Analysis 4 20 

19. TAD I/O pressure profiles for BFC computation: Analysis 4 20 

iv 



LIST OF ILLUSTRATIONS (Concluded) 


Figure Title Page 

20. TAD I/O turbulent kinetic energy profiles for BFC 

computation: Analysis 4 21 

21. Developing velocity profiles in the straight 2 x D exit 

section of Analysis 4 21 

22. Comparison of TAD outlet velocity profiles 23 

23. Comparison of TAD inlet pressure profiles 23 


v 


NOMENCLATURE 


Symbol 

C V>0 

Cm 


w 


D 

D 

E 

F 

k 

a 

i 

I 

N 

P 

P 


i] 


P.. 

U 


r. 

1 


R c 

S 
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e 


W u r ’ u 0 
U 


Definition 

constant (=0.09) in traditional k-e model turbulent viscosity 
expression 

spatially variant function in modified k-e model turbulent viscosity 

2 

expression (y t = C^pk /e) 

3/4 

constant in wall function correction (=k/C^q ), see equation (13) 

duct width 

dissipation term in the Reynolds stress transport equations 

east side cell, see Figure 12 

cubic polynomial in see equation (12) 

turbulent kinetic energy 

turbulence length scale 

turbulent length scale for flow over a plane 
north side cell, see Figure 12 
time averaged pressure 
turbulent energy production rate 
exit plane pressure assumed uniform 

production term in the Reynolds stress transport equations 

radial direction coordinate 

duct inner convex wall radius of curvature 

duct outer concave wall radius of curvature 

duct centerline radius of curvature (=[r.+r ]/ 2) 

1 o 

south side cell, see Figure 12 

momentum source terms, see equations (5) and (6) 

Reynolds stress components in cylindrical coordinates 
time averaged stream wise velocity 


vi 



Symbol 

U 

U r’ U 0 
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y eff 

n.. 

ij 
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e 

BFC 

FDP 

SSME 

TAD 


Definition 

cross- stream averaged stream wise velocity 

time averaged velocity components in cylindrical coordinates 

west side cell, see Figure 12 

distance from the duct inner convex surface measured in the 
radial direction 

empirical constant in curvature modified length scale, see equation 

( 1 ) 

energy dissipation rate 

dissipation term in the Reynolds stress transport equations 

von Karman constant 

absolute viscosity 

turbulent viscosity 

effective viscosity (=u t +p) 

pressure strain term in the Reynolds stress transport equations 
mass density 

angular position measured from the duct entrance plane 

body-fitted coordinates 

fully developed pipe or duct 

Space Shuttle Main Engine 

turn around duct 
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CONTRACTOR REPORT 


TURBULENT FLOW FIELD PREDICTIONS IN SHARPLY 
CURVED TURN AROUND DUCTS 

INTRODUCTION 


Performance demands of modern machinery require designs exhibiting both 
dimensional and operational economy. In applications which require the motion of a 
working fluid, these economic considerations have fueled a proliferation of devices 
containing complex flow passages through which fluids travel at extreme rates. 
Examples include flows in turbomachinery, nozzles, diffusers, and curved ducts. 

Of primary interest in this investigation are flows in sharply curved turn around 
ducts (TAD's) common in rocket engine design, including current and projected 
versions of the Space Shuttle Main Engine (SSME). 

Despite the design impetus provided by increased performance requirements, 
a detailed understanding of turbulent fluid motion with significant streamline curva- 
ture is inhibited by formidable obstacles to both experimental measurement and com- 
putational study. Since the system of motion governing equations for turbulent flows 
is essentially incomplete, advances in the computational prediction of turbulent field 
characteristics are tied to turbulence models requiring a sound empirical basis. In 
the case of sharply curved passages and extreme flow rates prevalent in modern 
equipment, the experimental basis required for detailed flow prediction is difficult to 
obtain with confidence. 

An example of a high rate, strongly curved flow field with little experimented 
basis on which to develop adequate turbulence models occurs in the SSME fuel side 
turbopump TAD. A simplified two-dimensional version of this TAD configuration is 
depicted in Figure 1. Although several studies of turbulent curved duct flow have 
been reported [1-7], none deals with curvatures or rates of the magnitude occurring 
in the SSME. The state of turbulence modeling in general is such that little confi- 
dence can be placed in extrapolation of the results of these investigations to a 
broader class of curved flows. 

The sensitivity of turbulent flow characteristics to even mild streamline curva- 
tures has been discussed by Bradshaw [8]. In this extensive review, a curvature 
modified length scale of the form 


t V r \ 
*= ( li60 _) to 


was proposed. In this relation g is an empirical constant of order 10, £ o is a length 

scale appropriate for plane flow, and subscripts following a comma indicate differen- 
tiation with respect to the specified variable (e.g., U 0 = 9U 0 /8r). Inclusion of 

this curvature modified scale into a general mixing length model of turbulent motion 
is hampered by the flow dependence of both g and l . Lack of extensive applicability 

is characteristic of all mixing length modeling strategies. 
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Launder et al. [4] have proposed a two-equation closure model with the length 
scale dependence obtained from solution of the equations for turbulent kinetic energy 
and the energy dissipation rate. The direct effect of curvature is introduced as a 
single term in the dissipation equation which is proportional to the turbulent Richard- 
son number. The required proportionality constant appears as a new empiricism. 

Recently, an extended k-e turbulence model, which uses algebraic approxima- 
tions for the Reynolds stress transport equations in order to develop an explicit 

9 

curvature dependent expression for in = C^k /e, has been reported by 

Pourahmadi and Humphrey [9,10] and applied to curved duct flows. The new expres- 
sion for accounts simultaneously for effects of streamline curvature and pressure 

strain with wall dampening. The development is a formidable task, calling upon a 
variety of traditional turbulence modeling approximations [11-15], The approach 
taken by Pourahmadi and Humphrey presents a promising explicit curvature model with 
roots, albeit approximate, in the Reynolds stress transport equations. Unfortunately, 
previous results based on this computational procedure appear to be based on a 
rather fundamental algebraic error which will be addressed in subsequent sections of 
this report. 

Another approach which has been pursued is to include curvature effects in a 
modified wall law [3,16]. This method is well suited for near wall treatment in coarse 
grid numerical calculations of turbulent boundary layers under adverse pressure 
gradients. Incorporation into computational procedures involving strongly curved, 
fully elliptic flow fields is certainly plausible. 

Because of the scarcity of reliable experimental data involving high rate flows 
in sharply curved passages, and the difficulty in extrapolating turbulent model pre- 
dictions, none of the above described curvature modified turbulence models can be 
objectively recommended. Each, however, presents an opportunity for computational 
investigation, refinement, and comparison with results of standard turbulence models 
as well as experimentation. 


OBJECTIVES 


The objectives of this investigation have a dual character. On one hand, com- 
putational predictions of turbulent flow properties in sharply curved TAD's are sought 
in order to provide additional information on the operational characteristics of devices 
incorporating such flows. Of specific interest are TAD flows in the SSME fuel side 
turbopump. As such, results of this investigation fall into the category of compu- 
tational analysis. A second view is more appropriate, however, with results of the 
investigation taken as presenting contributions to the base of computational experi- 
mentation needed to improve general models of turbulent flow with significant stream- 
line curvature. The data base building aspect of this study is considered to be the 
primary goal. 

In order to achieve this primary goal, three specific objectives were established. 

1) To estimate the general characteristics of turbulent flows in sharply curved 
TAD's using a standard k-e turbulence model on an abbreviated geometry as displayed 
in Figure 2. 
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2) To implement and test a specific turbulence model which explicitly incor- 
porates the effects of streamline curvature and pressure strain. 

3) To examine the elliptic character of the governing motion equations by 
implementing an orthogonal, body-fitted coordinate reference frame analysis on an 
extended duct geometry shown in Figure 3. 

Background information concerning each of the three stated objectives is provided in 
the following section of this report. 



e c 



Figure 3. Geometry and grid pattern for BFC TAD analyses: 20 x 80 grid. 
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BACKGROUND INFORMATION 


Equations governing the motion of steady, two-dimensional, turbulent flow of an 
incompressible, Newtonian fluid in cylindrical coordinates are given below [10]. 

Continuity : 


U r,r + U e,8 /r + V r=0 


( 2 ) 


r-Momentum: 


p[U r U r>r+ U 9 U r>6 /r - U 6 M = -P >r 4 <P e ff rU r,r>,r /r + <Peff U r, 9 /r) , 0 /r 


-P ef f« r /r2 - 2 Peff U 9,9 /r2 + S r 


(3) 


0-Momentum: 


P[U r U 9,r + U 6 U 9,6 /r + u r u e« = ‘V + <>‘eff rU 9,rV r + <“eff U 9, 9 /r) . 8 /r 


' 1, eff U e /r2 + 2ll eff U r,e /r2 + S 0 


(4) 


In equations (3) and (4) above, the Boussinesq assumption has been employed 
to model the Reynolds stresses in terms of mean velocity gradients and a turbulent 
viscosity component. The momentum source terms S r and S Q are given by the rela- 
tions below . 


s r = [u t r«yr) tr ] j9 /r + <V D r.r> ,r /r ‘ W‘ 


(5) 


S 6 = [p ( (2U r /r + U 6i6 /r)] >e /r 4 [p t (U r>6 - U e )/r] >r 4 „ t (U e _ p - U e /r)/r (6) 


In order to provide and solve for spatial variation of the turbulent viscosity 
y^, the k-e turbulence model as described by Launder and Spalding [17] was 

employed. This two equation model has been compared favorably to a variety of tur- 
bulence modeling schemes [ 18] . It has been tested on a wide range of flow con- 
figurations, exhibiting some degree of universality with an acceptable degree of com- 
plexity. Following this approach, y^ is determined by the relation 

y t = C pk 2 / e , (7) 
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where the turbulent kinetic energy and dissipation rate are determined upon solving 
the differential transport equations given below [10,19], 

Turbulent kinetic energy k : 


p[Uk 


u e k 


) /r] = ( ^eff rk ,r /a k\r /r + ( ^< 


sff 1 ^, e /o k ) : 


,/r 2 - 


P e + P, 


( 8 ) 


Energy dissipation rate e ; 


P [U r e ,r +U 9 e ,e /r] = <» 1 «ff re .r /a p ) .i. /r + ( ^ff e V r + C„ ie P^/k 


l eff ,r' e',r 


eff , 0 e' , 6 


el d ; 


- C e2 pe /k 


(9) 


The production term P d in equations (8) and (9) is given by the relation 

P d = “ t t 2 [‘ U r,r )2 + <U e,e /r)2 - U 9 (U r,0 /r + U 6,r> /r + U r< U r + 2U 9,e> /r2 


+ U r,6 U 6,r /r ) + <V r > 2 + (U 6,r> 2 + (U r,e /r)2) 


( 10 ) 


The numerical procedure used to solve the discretized form of equations ( 2) , 

(3), (4), (8) and (9) was provided by the proprietary PHOENICS computer code 
supplied and supported by CHAM of North America. A discussion of the philosophical 
and mathematical basis of this code is provided in Reference 19. A good presentation 
of the numerical procedures associated with a large family of CFD codes including the 
PHOENICS code is available in the text by Patanker [20]. 

The standard k-e turbulence model, the basis of which is presented above for 
cylindrical coordinate systems, was employed in computational analyses of simple 
geometry (Fig. 2) TAD flows. A 20 x 30 grid was employed in these analyses as 
displayed between planes BB and CC of Figure 3. In all simple TAD configuration 
analyses, inlet plane velocity, turbulent kinetic energy, and dissipation rate condi- 
tions were specified. A constant exit plane (CC in Fig. 3) pressure was prescribed 
and boundary condition closure was achieved by assuming zero velocity field gradients 
in the exit plane. 

The Pourahmadi and Humphrey [9,10] extension of the standard model was also 
implemented. As discussed above, this modified turbulence model incorporates the 
effects of streamline curvature and pressure strain in the fluid by using a variety of 
approximations to reduce the Reynolds stress transport equations to a set of algebraic 
equations. Because of the complexity of the algebraic manipulations employed, only 
a brief outline of the model derivation and implementation is presented below. The 
reader is referred to the original work in Reference 9 for a detailed development. 


The derivation starts from a high Renolds number form of the u.u. transport 
equation [9] characterized by the general form shown itr equation (11). ^ 
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( 11 ) 


D u.u. 

= P.. + 

Dt l] 


e.. + n.. + D.. 
i] i] i] 


The development then proceeds in the following stepwise manner. 

1) The u.Uj transport equations are reduced to a system of algebraic equations 
by invoking a variety of approximations including: 

A) Rotta's [12] return to isotropy principle inn. 


i] 


B) Launder, Reece and Rodi's [13] expression for the contribution to n. 
from interaction between mean strain and fluctuating velocities 

C) Daly and Harlow's [14] correction to n~ due to wall effects 


i] 


D) 

E) 


Gibson and Launder's [15] correction to n.. due to wall effects 
Rodi's [11] assumption that u.u./k is constant throughout the flow 

/»■ i -• 1 1 


field . 


2) The resulting system of algebraic equations is solved for the turbulent 
stresses. 

3) The Boussinesq approximation is then employed, transforming the solution 
into an algebraic relation for C^. 

4) The turbulent energy production term is rearranged into the form 


(k/ E )U e r = g(P d /c, 0 e /r, Uj.j )/<e; 


where g is a complex function of the quantities indicated, including the various par- 
tials of the mean velocity, U. .. 

*9 J 

5) The relation for (k/e)U from step 4 is substituted into the algebraic 

0 j r 

relation for from step 3. Simplification of the resulting expression leads to the 
cubic polynomial relation 

F(vTT~) = C 3/2 + a.C + a 9 C 1/2 + a„ = 0 (12) 

y y ly 2y o 


where a^a^a.^ = a 1 ,a 2 ,a 3 (P d /e, U Q /r, f, U.^) . 

6) The wall correction function f presented in Reference 13 and appearing in 
step 5 is modified to the form 


f = 


.3/2 


C w e 


I + (y/P) m 
y D-y 


(13) 
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thereby permitting asymmetric wall corrections to the pressure strain. In this 
expression, the asymmetry measure m appears as a new empirical constant. 

7) Following Bradshaw [8], the expression for dissipation at near wall nodes 

is modified by a multiplicative factor of 1/[1±BU ./(U Q r)] in order to introduce the 

o u,r 

influence of extra strain curvature on the magnitude of the near wall turbulence 
length scale. This expression also contains a new empirical constant, 3. 

Suggested values for constants m and 3 are 7.95 and the order of 10 respectively. 

Although the complex development terminating at equation (12) is quite impres- 
sive, an unfortunate algebraic error is made in establishing the roots of this equa- 
tion. Instead of a single positive real root for as predicted in References 9 and 

10, multiple real roots, not necessarily positive, can (and did) occur in practice. 

The implemented version of the Pourahmadi and Humphrey modification included a 
corrected method of solving for the curvature dependent values of t/C^ . Adding to 

the difficulties associated with this implementation was the need to prescribe a selec- 
tion criterion for regions in which multiple roots occurred. This problem was never 
resolved with complete success as will be discussed more fully in the next section of 
this report. 

In order to more fully understand the elliptic character of developing TAD 
flows, analyses utilizing an orthogonal body-fitted coordinate (BFC) reference system 
on an expanded duct geometry were performed. Two such analyses were performed; 
one on the full geometry, 20 x 80 grid system depicted in Figure 3, and a second on 
a reduced geometry 20 x 60 grid depicted between planes AA and CC of the same 
figure. In each case the 20 x 30 grid configuration of the simple TAD analyses was 
preserved in the region bounded by surfaces BB and CC in Figure 3. In this manner 
comparison of results using both simple TAD configurations and expanded BFC geo- 
metries was facilitated. 

Numerical computations employing the grid configuration of Figure 3 were per- 
formed utilizing the BFC option of the PHOENICS computer code described in Refer- 
ence 21. Although potentially more accurate non-orthogonal grid specifications could 
have been implemented within the framework of the PHOENICS code, direct comparison 
with results of simple cylindrical system analyses would have then become more diffi- 
cult. An excellent comprehensive review of methods to numerically generate curvi- 
linear coordinate systems with coordinate lines coincident with boundary segments is 
given in Reference 22. 


RESULTS 


Results of this investigation can be separated into three categories dependent 
upon coordinate reference frame and turbulence model employed. Analyses 1, 2, and 
6 described in Table 1 were performed on simple TAD geometries using a cylindrical 
coordinate system and a standard k- e turbulence model [17]. In analysis 3 the 
modified k-e model of Pourahmadi and Humphrey [9,10], with the correction discussed 
in the previous section, was employed. The simple TAD geometry displayed in 
Figure 2 was maintained. Analyses 4 and 5 employed orthogonal BFC reference 
systems with a standard k-e turbulence model. The grid system depicted in Figure 3 
was used in analysis 4. The same grid pattern, restricted to the corresponding 
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analysis geometry, was employed in all cases. In addition, the flow Reynolds number 
6 

was fixed at 10 and a curvature ratio R c /D = 1 was prescribed. These values 

obviously refer to a high rate flow with extremely sharp curvature. A complete sum- 
mary of the computational analyses of this study is presented in Table 1. 

TABLE 1. ANALYSES SUMMARY 


Analysis 

Reference 

Figure 2 

Grid 

Inlet U 

Turbulence 

No. 

Frame 

I/O Planes 

Profile 

Model 

1 

Cyl. 

BB/CC 

20 x 30 

Plug 

k — e 

2 

Cyl . 

BB/CC 

20 x 30 

FDP 

k-e 

3 

Cyl . 

BB/CC 

20 x 30 

Plug 

Modified 






k-e 

4 

BFC 

AA/EE 

20x80 

Plug 

k-e 

5 

BFC 

AA/CC 

20x60 

Plug 

k-e 

6 

Cyl. 

BB/CC 

20 x 30 

Plane BB 

k-e 





from 5 


Analyses Constants: 

Re = 10 1 

5 R /D 

c 

= 1.0 


Turbulence Inlet Conditions 
For Plug Inlet U Profile : 

k/U 2 = 0.005 





z/ (U 3 /D) = 

(0.005) 3/ 2 

/0.01 

Exit Plane 

Pressure: 

Uniform 

across exit 

plane 



In Figures 4, 5, and 6 are displayed developing profiles of velocity, pressure, 
and turbulent kinetic energy respectively, for the case of a simple TAD geometry 
with inlet plug flow. In response to the strong cross stream pressure gradients 
evident throughout the flow domain, the velocity profiles are observed to become 
increasingly asymmetric. The most substantial pressure adjustments occur within the 
final 30 deg of duct travel, and in response there is a significant acceleration of flow 
near the outer wall with a corresponding deceleration near the inner, convex, duct 
wall. The net effect is a nearly linear velocity profile very near the exit plane 
reminiscent of a solid body rotation. Quantitatively, the near exit plane velocity 
varies from just under 25 percent to somewhat over 140 percent of the average flow 
velocity. In the entrance region of the duct, the plug flow characteristic is observed 
to persist virtually unchanged past the midplane ( 0 = 90 deg) of the duct. The large 
cross stream pressure gradient predicted for the TAD inlet region together with the 
specified plug flow entrance velocity profile is suggestive of a very large swirl 
chamber feed to the TAD . This type of configuration is of course not prevalent in 
rocket engine or turbomachinery applications. 

The turbulent kinetic energy profiles displayed in Figure 6 exhibit rapid decay 
in the plug flow core region as expected. An initial decay in the near wall turbulent 
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kinetic energy is followed by a significant rise in the accelerating concave wall region, 
and a continuing gradual decay in the decelerating convex wall region. Near the exit 
plane, the turbulent energy profile develops the knee observed [1,2] and predicted 
[10] in previous curved duct investigations conducted at smaller values of Re and R c /D . 



Figure 4. Streamwise velocity profiles for simple TAD geometry 
with plug flow inlet: Analysis 1. 



<pU 2 /2) 

Figure 5. Pressure profiles for simple TAD geometry with 
plug flow inlet: Analysis 1. 
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Results of a simple TAD geometry analysis with a fully developed turbulent duct 
flow inlet condition are displayed in Figures 7, 8, and 9. The general shape charac- 
teristics and scale of the pressure profiles in Figure 8 are similar to the previous 
plug flow inlet case. The development sequence for the velocity profiles is also simi- 
lar to the analysis 1 case with little profile adjustment prior to the 90 deg midplane 
of the TAD , and significant velocity gradient increases and profile linearization in the 
rapid adjustment exit region. The salient feature displayed in Figure 7 is the appear- 
ance of a recirculation region near the inner wall of the duct. This region extends 
back somewhat over 6 deg from the exit plane and occupies approximately 4 percent 
of the duct cross section at the exit plane. The occurrence of a negative velocity 
region in the exit plane presents certain difficulties involving computational error as 
discussed in Reference 21. Despite this problem, flow separation from the inner con- 
vex wall is predicted to occur between 0 = 171 deg and 0 = 174 deg. This phenomena 
is not observed to occur in less sharply curved channels with smaller values of the 
flow Reynolds number [1,2,6]. 

The turbulent kinetic energy profiles depicted in Figure 9 differ substantially 
from the plug flow inlet case both in shape and development history. The initial 
profile shape is symmetric and of the same magnitude as displayed for the analysis 1 
entrance profile. Because of the entrance cross stream velocity gradients present 
in the FDP entrance case, however, the turbulent kinetic energy is predicted to 
increase significantly in the flow direction away from the convex wall . This is in 
sharp contrast to the rapid post entrance decay of k exhibited in Figure 6 for the 
plug flow inlet case. Near the exit plane, cross stream velocity gradients and the 
energy dissipation rate away from the convex wall increase rapidly in the flow direc- 
tion, and the increase in turbulent energy is arrested. In the vicinity of the convex 
wall there is a substantial adverse pressure gradient leading to deceleration and, as 
the exit plane is approached, flow separation. This is accompanied by a sharp 
decrease in the inner wall region values of k. 

Direct comparison of the downstream velocity and pressure profiles predicted in 
analyses 1 and 2 are displayed in Figures 10 and 11. The similarities in each of the 
profile pairs is obvious. The near exit plane pressure adjustment is predicted to 
occur more suddenly in the case of plug flow inlet conditions. It is evident that the 
inner wall characteristics of the near exit plane flow are substantially affected by 
changes in the inflow boundary conditions. 

In the third analysis of the computational sequence, the modified k-e model of 
Pourahmadi and Humphrey [9,10] was implemented in the user definition section of 
the PHOENICS routine. Severe difficulties were initially encountered in attempts to 
obtain reasonable solutions. Further investigation revealed the equation (12) rooting 
error described previously. 

In the context of the PHOENICS implementation, the Pourahmadi and Humphrey 
[9,10] modification for a specific finite difference cell P requires manipulation of the 
staggered grid velocity vectors indicated in Figure 12. In addition, grid distances 
from each of the indicated neighboring cells to the P cell boundaries are needed 
together with estimates of k and e at P. The entire manipulation was computationally 
intense, often tripling the clock time of analysis. 

Despite the increased computational effort, convergence difficulties were 
repeatedly encountered. The source of the problem was the behavior of the poly- 
nomial F(v^r~ ), the roots of which indicate candidate values of C . Typical F(/C ) 

V PM 

loci are depicted in Figure 13. These polynomial traces were arrived at by 


12 



0.0 0.5 1.0 1.5 


P~Pe 

(p 02/2) 

Figure 8. Pressure profiles for simple TAD geometry with FDP 

flow inlet: Analysis 2. 
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Figure 9. Turbulent kinetic energy profiles for simple TAD 
geometry with FDP flow inlet: Analysis 2. 













computationally testing the modified k-e turbulence model on a simple TAD geometry, 
plug inlet flow problem. Analysis 3 as described in Table 1 was initiated at the 
standard k-e model solution value of analysis 1. An intermediate flow field result 
was determined by arbitrarily selecting the principal theta transformation root of the 
cubic polynomial F(v / 5^) whenever multiple roots occurred. The trigonometric or 

theta transformation method for determining roots of a cubic polynomial is described 
in most mathematical handbooks (see for example Reference 23). 

It is evident from Figure 13 that for 6 = 168 deg, the near concave wall 
behavior of F was that of a cubic polynomial with multiple real roots. This behavior 
was not isolated, occurring irregularly throughout the flow field and frequently in 
near wall regions. Unfortunately it is in these near wall regions that the modified 
k-e model was expected to yield its most meaningful results [9,10]. The behavior 
of the F(vtr^) polynomial at y/D = 0.9875 (the node closest to the concave wall) was 

particularly disturbing. As indicated in Figure 13, the trace of F was virtually 

colinear with the C axis for -0.2 < C <0.2. Three distinct real roots were deter- 

y y 

mined at this nodal position, however, the ill natured character of F is apparent, 
making the computational process extremely susceptible to even small numerical errors. 


Near exit plane values of were determined as shown in Figure 14. It should 
be noted that the line of = 0.09, the standard k-e model value, occurring at 0 = 

150 deg was an assigned default value selected whenever the cross stream velocity 

gradient became small. Small values of U„ led to extreme variation of the Ffv'C ) 

0 , r y 

polynomial, hence the default selection in these regions. Pourahmadi and Humphrey 
[9,10] also note the ill behavior of in the core region of the duct. They attribute 

this behavior to a breakdown in the assumption that uTuT/k is constant in this region. 

Their reasoning is less than compelling. It is entirely possible that the problem of 
determining reasonable values for C is completely numerical, and that consideration 


of a more primitive implicit form of the C relation might yield stable results across 
the entire duct cross section. y 


Recognizing that the analysis 3 results are only intermediate, lacking a com- 
pletely acceptable method of selecting C^, it is useful and interesting to observe 

position dependent characteristics of the variant . In Figure 14, for 0 = 174 deg, 

is observed to take on values substantially less than the standard value of 0.09 

in near wall regions, rising to nearly 0.09 just outside the convex wall region and 
then falling to near zero in the flow core, rising again to near the standard value and 
well beyond just outside the concave wall region. Figure 15 displays the near concave 
wall behavior of throughout the flow domain. It is evident that the near wall 

node value of remains essentially constant at around 0.03 throughout while the 

adjacent node value of took on default values near the TAD entrance, rising 

dramatically to a value of nearly 0.4, and then falling steadily to a value in excess 
of 0.1 before rising dramatically near the duct exit. The turbulent kinetic energy 
profiles displayed in Figure 16 exhibit the same qualitative character as the unmodi- 
fied k-e model profiles of Figure 6 except that near outer wall k values are greatly 
magnified and exhibit a knee as the wall is approached. Details of the near convex 
wall comparisons of turbulent kinetic energy are shown in Figure 17. Both the 
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Figure 16. Turbulent kinetic energy profiles for simple TAD 
geometry using the modified k-e model: Analysis 3. 
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Figure 17. Comparison of near wall turbulent kinetic energy predictions 



standard and modified k-e models predict the same k values at the adjacent wall node, 
however, modified model kinetic energy predictions were somewhat smaller than 
standard at the neighboring y node near the duct exit. Estimated velocity and pres- 
sure profiles obtained using the modified k-e model were virtually identical to stan- 
dard model results. 

Two separate BFC computational analyses were performed in order to explore 
the elliptic character of the TAD flow field. These are described in Table 1, analyses 
4 and 5. Before describing results of these studies, several general observations are 
warranted. First, the early version of the SSME TAD geometry depicted in Figure 1 
suggests a computational model with a short entrance region and a relatively long 
straight exit region. In Figure 3 this would indicate flow entrance at plane EE and 
exit at plane AA. Several attempts at such an analysis using the BFC option of the 
PHOENICS computer code were unsuccessful. This does not necessarily indicate a 
shortcoming in the PHOENICS /BFC solution technique. It is likely that more extensive 
experience with judicious selection of the relaxation factors and other heuristic tuning 
parameters inherent in the PHOENICS code, and indeed any CFD code with like 
capabilities, would have led to solution convergence with the flow orientation 
described above. Because of the limited time available to gain computational experi- 
ence, and because it was found that extension of the straight duct inlet zone greatly 
enhanced convergence characteristics, the suggested SSME TAD flow orientation was 
reversed, with the inlet prescribed in plane AA and exit at plane EE. 

A second observation concerning the BFC computations involves the degree of 
solution convergence at analysis termination. Using the absolute residuals sum cri- 
terion described more fully in the PHOENICS documentation, the convergence rate at 
termination had slowed to approximately 0.1 percent of total absolute residuals sum 
per solution sweep. This rate was achieved after approximately 30 clock hours of 
computation on the Perkin-Elmer 3220 computer system utilized during this study. 
Although clock run time is a poor absolute measure of computational efficiency, rela- 
tive to simple TAD configurations, the BFC analyses required an order of magnitude 
more clock run time and at best achieved an absolute residuals sum convergence 
measure a full order of magnitude larger than the simple TAD analyses. Noting that 
the BFC grid patterns were approximately three times the size of the simple cylin- 
drical system TAD grids, it is obvious that the BFC analyses were relatively high 
overhead computations. In addition, because calculation was interrupted before 
achieving a proportionately small convergence measure, results of the BFC analyses 
are considered to yield only qualitative predictions of flow field characteristics. 

Figures 18, 19, and 20 display results of the analysis 4 BFC computations for 
TAD inlet /outlet velocity, pressure, and turbulent kinetic energy profiles respectively. 
These calculations were performed using the full grid pattern displayed in Figure 3 
with flow inlet at plane AA . It is evident in Figures 18 and 19 that BFC calculated 
TAD inlet flow field features are radically different from the inlet boundary assump- 
tions of the simple duct analyses. The inlet velocity profile resembles the 1/r pro- 
portionality of ideal fluid flow except in the near wall regions while the outlet profile 
appears very much like a fully developed turbulent duct flow. Obviously the anti- 
symmetric character of the predicted TAD outlet pressure profile is quite different 
from the assumed uniform exit plane distribution of the simple geometry analyses. 
Because of the centrifugally induced cross- stream pressure gradient at the TAD exit, 
the duct like flow profile at plane CC developed into the stretched, r proportional 
profile reminiscent of rigid body rotation in the 2 x D straight exit duct. Predicted 
developing velocity profiles in the exit duct are displayed in Figure 21. In the exit 
section, the pressure adjusts to the assumed uniform condition at exit plane EE. The 
core of the TAD exit plane k profile shown in Figure 20 resembles the core suppressed 
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profile of the simple geometry analysis with plug flow inlet conditions displayed in 
Figure 6. No inner wall knee is observed in the BFC k predictions, however, outer 
wall turbulent kinetic energies are reduced in the BFC analysis. 

In order to reinforce the extended geometry predictions of analysis 4, a com- 
parison test was devised. BFC analysis 5 was run on the reduced configuration 
between planes AA and CC . Plug flow inlet conditions were prescribed at plane AA 
and uniform exit plane pressure was assumed. Results of this computational proce- 
dure at the TAD entrance plane BB were used as boundary conditions for a simple 
geometry analysis using a cylindrical system with the same uniform exit plane pressure 
assumed. Predictions of analysis 5 and 6 should have been identical, at least theo- 
retically. Unfortunately, there was a substantial discrepancy in the predicted flow 
field characteristics as demonstrated in Figures 22 and 23. Even considering the 
qualitative nature of the BFC analysis, the comparisons were disturbing. Cylindrical 
system predictions of cross- stream velocity gradients at the exit plane were generally 
much larger than results using BFC analysis as is evident in Figure 22. In addition, 
a rather large recirculation zone was predicted using the traditional cylindrical coor- 
dinate system analysis. Inner wall flow separation was predicted approximately 24 deg 
from the exit plane, with the counterflow velocity zone extending into the fluid over 
7.5 percent of the duct width at the exit plane. By contrast, the near inner wall 
node velocity at the exit plane, obtained using BFC computations was over 40 percent 
of the average velocity. Pressure profiles from the two analyses at the TAD entrance 
plane were not even similar in character, as displayed in Figure 23. 

If the standard cylindrical system analysis is presumed correct, based on an 
extensive history of computational success and better convergence characteristics, 
then the reliability of the current PHOENICS BFC option must be suspect. This con- 
clusion must, however, be tempered by the fact that acceptable convergence was 
attained with cylindrical system analysis while only run time constrained convergence 
was achieved using BFC's. In addition, in order to prevent divergence of the 
PHOENICS based computational procedure involving BFC's, it was necessary to neglect 
certain curvature contributions in the finite difference approximations. The rationale 
for neglect of these terms is more fully described in the PHOENICS BFC documenta- 
tion [21], however, for curvature dominated flows, such as that occurring in the 
SSME TAD, neglect of such contributions is dissatisfying. Only by further computa- 
tional study and experimental comparison can deficiencies in the approximation 
procedure be adequately quantified. 


CONCLUSIONS 


1) The curvature modified k-e turbulence model developed by Pourahmadi and 
Humphrey [9,10] cannot be utilized as presented due to algebraic errors in the 
F(vtT) rooting scheme. 

2) A corrected version of the Pourahmadi and Humphrey [9,10] curvature 
modified turb ulen ce model has potential if an adequate method of selecting the multi- 
valued root can be determined. 

P 

3) Neglect of curvature related terms within the BFC option of the general 
purpose CFD code PHOENICS yields results inconsistent with the predictions of 
standard orthogonal systems analyses for sharply curved internal flows. 
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Figure 23. Comparison of TAD inlet pressure profiles 
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